Self— Similar Perturbation Theory 

V. I. Yukalov 1,2 and E. P. Yukalova 2,3 

1 Bogolubov Laboratory of Theoretical Physics 
Joint Institute for Nuclear Research, Dubna 141980, Russia 

^ ' 2 Instituto de Fisica de Sao Carlos, Universidade de Sao Paulo 

Caixa Postal 369, Sao Carlos, Sao Paulo 13560-970, Brazil 



OS 



q \ ^Department of Computational Physics 



in 



- 1—1 

X 



Laboratory of Computing Techniques and Automation 
Joint Institute for Nuclear Research, Dubna 141980, Russia 



Abstract 



> 

in 

CO 

y—i 

0\ 

A method is suggested for treating those complicated physical problems for which 
exact solutions are not known but a few approximation terms of a calculational 
Ph. algorithm can be derived. The method permits one to answer the following rather 



delicate questions: What can be said about the convergence of the calculational 
procedure when only a few its terms are available and how to decide which of the 



initial approximations of the perturbative algorithm is better, when several such 
initial approximations are possible? Definite answers to these important questions 
become possible by employing the self-similar perturbation theory. The novelty of 
this paper is in developing the stability analysis based on the method of multipliers 
and in illustrating the efficiency of this analysis by different quantum-mechanical 
problems. 
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1 Introduction 



Realistic physical problems can practically never be solved exactly quantum field theory 
being among these the most difficult problems. Then one has to resort to some approxi- 
mations. One such a reliable approximation is the Gaussian-effective-potential approach 
[1-5]. This approach contains the leading 1/N result as its formal N — > oo limit and 
also the one-loop result as its formal H — > limit. Thus, it encompasses both of the 
other popular approaches to effective potentials, transforming their leading-order expan- 
sions and displaying a much richer structure [5,6]. At zero temperature, the Gaussian 
effective potential can be described as a variational approximation to the vacuum energy 
density constructed with the Gaussian trial wave functions shifted by a constant classical 
background field [7]. The finite-temperature Gaussian effective potential is just the ap- 
proximate free energy in the self-consistent harmonic approximation [7,8]. In this way, the 
Gaussian effective potential in quantum field theory is a direct analog of the variational 
ground-state energy in quantum mechanics, so that the quantum-mechanical language 
can be directly transcribed to the quantum-field-theory language [1]. The difference is 
that variational methods in quantum mechanics permit one to derive several nice com- 
parison theorems [9-13] for the ground-state energy, which is rather difficult, if possible, 
in quantum field theory. 

To obtain corrections to a chosen initial approximation, one has to develop some per- 
turbative algorithm about the given approximation. Thus post-Gaussian corrections to 
the Gaussian effective potential may be obtained [14,15], allowing the variational param- 
eter to change from one order to the next for the expansion to yield convergent results 
[14-17]. In the same spirit, one can calculate the effective potentials or energies as well 
as wave functions [18,19]. Different variants of such an approach are called in literature 
by various names, as modified perturbation theory, renormalized perturbation theory, 
optimized perturbation theory, controlled perturbation theory, oscillator-representation 
method, and so on [20-31]. One often dignifies this kind of methods as nonperturbative, 
implying that they result in modified, or resummed, expansions having a more complicated 
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structure than simple weak-coupling power series. However, the term "nonperturbative" 
is, to our mind, not only cumbersome but somewhat misleading. This is because from 
the mathematical point of view any regular procedure producing a sequence of approxi- 
mations is perturbation theory, irrespectively of what initial approximation is chosen and 
what additional conditions in the course of calculation are imposed. The role of such 
additional conditions is to reorganize, by defining control functions, a given series into 
another one with better convergence properties [20-22] , which is necessary because of the 
divergence of standard weak-coupling power series [32]. Note that the renormalization 
scale, appearing in perturbative expansions for observable quantities of field theories, can 
also be treated as a control function that may be defined by some additional conditions 
[27,28,33] or invoking information from experiment [33,34]. 

The convergence of modified perturbation theory can be treated in three ways. The 
straightforward case is when the model considered is simple enough to permit the eval- 
uation of perturbative terms of arbitrary large orders. This happens for zero- and one- 
dimensional anharmonic oscillators [35-37]. But clearly one cannot expect to have the 
same luck for less trivial models. 

Another possibility to check convergence is when the exact or very accurate numerical 
solution of the problem is available. This is so for many quantum-mechanical problems. 
Then one may directly compare the calculated perturbative terms with the known accurate 
solutions. If these terms approach the latter, one tells that the perturbation procedure is 
convergent. Since in this case one deals with a finite number of approximate terms, it is 
more correct to say that one observes numerical convergence or local convergence. It is 
this type of convergence that one deals with in any calculational procedure stopping at 
a finite-order approximation. Certainly, the luxury of possessing for comparison exact or 
almost exact numerical results is rather rare for realistic physical problems. 

The numerical or local convergence may be also observed when one is able to find 
quite a number of perturbative terms, say several tens or hundreds of them, so that 
the calculational procedure saturates at some value that changes less and less with the 
increasing order of approximation. Again, it is not too realistic to hope to reach such a 
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saturation with tens or hundreds of perturbative terms in complicated problems, such as 
field theory, usually allowing just a few terms. 

What then can be said about the reliability of the results obtained by means of a 
perturbation algorithm, even if this is an optimized or modified perturbation theory, 
when one has in hands only a few terms and no exact answers are known? It may 
happen that the subsequent perturbative results quantitatively are not close to each other, 
then which of them is preferable? For example, optimized perturbation theory has been 
systematically tested using different ansatze for a range of different physical quantities in 
quantum electrodynamics and quantum chromodynamics [33]. Although in some specific 
applications this theory has met with success, in the majority of tests the optimization 
procedure was not successful. Typically, though the sign of the coefficients are correctly 
predicted, the optimization ansatze give values that are an order of magnitude too large 
or too small. Thus one comes to the conclusion that it is difficult to be optimistic about 
the general usefulness of optimized perturbation theory [33]. Even the more dramatic 
situation is when the neighbour-order results differ qualitatively, thus yielding different 
physics. For instance, for (2 + l)-dimensional scalar theory, the analysis of the Gaussian 
effective potential shows a first-order phase transition, in contrast with the post-Gaussian 
approximation indicating a phase transition of second order [15]. Another example is the 
so-called autonomous version of the scalar field theory, which exists in the Gaussian 
approximation but does not survive, becoming unstable, in the higher orders of the post- 
Gaussian expansion [15]. 

Thus, for any calculational algorithm, including all variants of modified or optimized 
perturbation theory, there exists the general question: "How can one trust to a pertur- 
bative procedure when only a few first approximations are available and neither exact 
solutions, nor accurate numerical calculations, nor detailed experiments are known for 
comparison?" To say this in other words: "Is it possible to formulate a calculational algo- 
rithm that would be equipped by an internal self-consistent criterion allowing to estimate 
the validity of the obtained approximations?" In the present paper we aim at suggesting 
a positive answer to this question. The novelty of this paper is twofold: First, we de- 
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velop the stability analysis, based on the method of multipliers, so that it really becomes 
possible to check the validity of the calculated approximations. Second, we illustrate 
the efficiency of the developed analysis by a number of quantum-mechanical problems. 
Among the latter, we consider the zero-dimensional (p 4 theory and anharmonic oscillators 
with integer powers. Although these models have been treated earlier by the self-similar 
approximation theory, however the stability analysis, as developed in this paper, have not 
been applied to them. The main part of the paper is devoted to the problems that have 
not yet been treated by the self-similar perturbation theory. These are the Hamiltonians 
with power-law potentials, having arbitrary noninteger powers, and a Hamiltonian with 
a logarithmic potential. In this way, the results presented in the paper are new. 

2 Basic Formulas 

Our consideration is based on the self-similar approximation theory [38-45]. This ap- 
proach is more general than modified perturbation theory because of the possibility to 
reformulate calculational procedure as the problem of analyzing the evolution equation 
for a controlled dynamical system, thus, permitting us to employ powerful techniques of 
dynamical theory and optimal control theory. We shall not repeat here the mathemati- 
cal foundations for the self-similar approximation theory, which have been expounded in 
detail in earlier papers [38-45]. We will only mention the main idea of the approach and 
delineate its scheme, necessary for the following illustrations; but we shall concentrate 
on the stability problem that has not yet been considered carefully, the problem which is 
pivotal for the motivation of this paper, and which would make it possible to answer the 
questions formulated in the Introduction. 

Assume that we are looking for a function f(x) of a variable x. For simplicity, 
we imply that the function and variable are real, though the whole procedure can be 
straightforwardly extended to the case of complex functions and variables. Let the sought 
function f(x) cannot be found exactly but only its perturbative expressions Pk(x), with 
k = 0,1,2, .. ., can be obtained in the asymptotic vicinity of some point, say x — 0. That 
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IS, 



f(x) ~ p k (x) 




0). 



(1) 



The standard form of Pk(x) is a series in powers of x, because of which x is called the 
expansion parameter. This can be, e.g., a coupling parameter. Such series are practically 
always divergent. If a number of terms Pk(x), about tens or more, were known, one 
could invoke some resummation techniques for ascribing a meaningful value to a divergent 
series [46]. But these techniques are useless when only a few perturbative terms are 
available. With the knowledge of just a few terms, one could find the limit of a sequence 
if a recurrence relation between subsequent terms were known. However, such recurrent 
relations are not merely difficult to discover but they have sense only for convergent 
sequences. Really, for an asymptotic series Pk{x), meaningful solely for asymptotically 
small x — > and diverging for any finite x, a recurrence relation between the terms 
Pk(x) and pk+i(x) cannot be defined since the limit of Pk(x), as k — > oo, does not exist. 
Nevertheless, if the sequence of approximations {pk(x)} is obtained from one given system 
of equations, corresponding to a physical problem under consideration, by means of the 
same perturbative algorithm, the set of found terms pk(x) does contain information about 
the sought function f(x). But this information, because of divergence of the sequence 
{pk{x)}, is hidden, or we can say that it is enciphered, encoded. 

The first thing we have to do in order to decipher, to decode the hidden information 
is to reorganize the sequence {pk(x)} to another sequence that would be convergent at 
finite x. In the process of this reorganization, or renormalization, that can be symbolically 
denoted as 



we introduce auxiliary functions u = Uk(x) so that the sequence {fk(x)} of the terms 



be convergent [20]. Because of the role of the functions Uk(x) in controlling convergence, 
they can be named the control functions. For a convergent sequence {fk(x)} a limit f*(x), 
as k — > oo, exists, which represents the sought function f(x). Now the problem of finding 
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U{p k (x)} = {F k (x,u)} , 



(2) 



fk(x) = F k (x,u k (x)) 



(3) 



a recurrence relation, or mapping, between subsequent terms fk(x) and fk +i (x) can be 
meaningful. 

The main idea in finding a relation between the terms of a sequence {fk(x)} is to 
try to construct a dynamical system with discrete time, being the approximant number 
fc = 0,l,2,..., so that the motion from fk(x) to fk+ n (x) would correspond to the evolution 
of this dynamical system. For this purpose, we need, first, to pass from the sequence 
ifk(x)} to a set of endomorphisms, which is done as follows. Define the expansion function 
Xk{<p) by the equation 

F Q (x,u k (x)) = (p , x = x k ((p). (4) 
Introduce an endomorphism y k by the transformation 

Vk(<p) = fk(x k ((p)) , (5) 

where, by definition (4), yo((f) = (p. As is evident, the sequences {ykif)} and {fk(x)} are 
bijective, and the limit y*((p) of the former sequence is in one-to-one correspondence with 
the limit f*(x) of the latter sequence. The limit y*(ip) of the sequence of endomorphisms 
{ykif)} is nothing but a fixed point for which 

y fc (y»)=y». (6) 

A family of endomorphisms in the vicinity of a fixed point satisfies the equation 

Vk+niv) = Vk(y n (<p)) , (7) 

which can be easily shown by checking that, for tp — > y*, Eq. (7) reduces to the identity 
y* = y*. In physical parlance, Eq. (7) is termed the property of self-similarity, the usage 
of which justifies the adjective self-similar, as is employed with respect to perturbation 
theory in the title of this paper. In the mathematical language, the equations y k + n = yk-y n 
and yo — 1 are the semigroup properties. A family of endomorphisms {y k \ k — 0, 1, 2, . . .}, 
equipped with the semigroup properties, forms a cascade, that is, a dynamical system with 
discrete time. Since, by construction, the cascade trajectory {yk(<p)} is bijective to the 
approximation sequence {fk(x)}, we call {y k \ k — 0, 1, 2, . . .} the approximation cascade. 
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To consider time evolution, it is convenient to deal with continuous time. For this, we 
may embed the cascade into a flow, 

{s/ fc |fc = 0,l,2,...}e{s/ t |f >0}, (8) 

which implies that the flow trajectory passes through all points of the cascade trajectory, 

yt(<p)=y k (<p) (t = fc = 0,l,2,...), (9) 
and the same semigroup property, as for the cascade, is true for the flow, 

yt+Av) = ytivM) • (10) 

The flow {y t \ t > 0} embedding the approximation cascade is called the approximation 
flow. 

From the group relation (10), the Lie equation 

d 

— y t ((p) = v(y t (ip)) (11) 

immediately follows, in which v(y) is a velocity field. We need to find a fixed point y*(ip) 
of the approximation flow whose motion is given by the evolution equation (11). Consider 
the latter starting from the point y k {tp) ait — k and moving to an approximate expression 
for the fixed point y k+1 (f, t) = 7/£ +T (<^) occurring at the time t — k + r. Then, integrating 
Eq. (11), we come to the evolution integral 



f U '^ = r, (12) 
Jyu v(y) 



in which y* k+1 = y* k+1 ((p, r) and y k = y k (tp). 

To make the evolution integral (12) useful, we have to concretize the flow velocity. The 
latter can be done by the Euler discretization giving for the velocity v(yk((p)) = Vk((p) the 
form 

d 

v k (ifi) = F k+1 (x,u k ) - F k (x,u k ) + (u k+1 - u k )—— F k (x,u k ) , (13) 

where x = x k (ip) and u k = u k (x k (ip)). The discretized velocity (13) may be termed the 
cascade velocity. To obtain an explicit expression for the latter, we need yet to define 
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the control functions Uk(x). This can be realized being based on the following reasoning. 
As is clear from the Lie evolution equation (11), at the fixed point y*, when y t (y*) = y*, 
the flow velocity is zero, v(y*) = 0. Hence, closer we are to a fixed point, smaller is the 
absolute value of velocity. This suggests that the fastest convergence to the fixed point is 
achieved if control functions provide the minimum of the velocity modulus 



My?) I < \F k+1 (x,u k ) - F k (x,u k )\ + 



d 

(u k+1 - u k )—^F k (x, u k ) 



(14) 



To minimize the second term in the right-hand side of the latter inequality, we may set 

d 

(u fc+ i - u k )-— F k (x,u k ) = . (15) 
ou k 

The principal usage of Eq. (15) is as follows. We are looking for a solution of the equation 

d 

—F k (x, u k ) = , U k+l ^ u k , 

which is labelled as the principal of minimal sensitivity [27], and which defines the control 
function u k {x). If the latter equation has no solution, we put u k+ i = u k . If this equation 
possesses several physically meaningful solutions, we must take that of them which min- 
imizes the cascade velocity (14). In the other way, we could set to zero the first term in 
the right-hand side of Eq. (14), thus coming to the principle of minimal difference [20]. 
But then we should invoke some additional condition for defining u k+1 . In what follows, 
we use the fixed-point condition (15), so that the cascade velocity (13) becomes 

v k ((p) = F k+1 (x, u k ) - F k (x, u k ) , (16) 

where x = x k (ip) and u k = u k (x k (ip)). 

In this way, taking into account Eqs. (3)-(5), we rewrite the evolution integral (12) 
in the form 

Jfk v k {(p) 

in which = f k (x, r) is a A;-order self-similar approximant for the sought function f(x) 
and f k = f k (x). The parameter r in Eq. (12) is an effective minimal time necessary for 
reaching a fixed point y k+1 - If no constraints are imposed on the properties of the sought 
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function, then the minimal time is the minimal number of steps, that is r = 1. When 
some additional conditions are imposed, which the sought function has to satisfy, then 
t = Tk(x) is to be treated as another type of control functions defined by these conditions. 
One more possibility of fixing r is connected with the problem of stability, as is described 
below. 

The problem of stability is extremely important. This is just what permits us to decide 
which of the found self-similar approximants are trustworthy and which of them are better 
then others. Let a set of mappings {y k (<p)} be given. Stability analysis is based on the 
concept of the local mapping multipliers 

d 

»k(<p) = Q^yk(v) ■ ( 18 ) 

When 

\M\ < 1 (19) 

for a fixed k, one says that the mapping at a fc-step is locally stable with respect to 
the initial point if. If \/i k (ip)\ = 1, one tells that the mapping is locally neutral. In our 
case, the sequence {y k (ip)} is the trajectory of the approximation cascade. Therefore, if 
inequality (19) holds true, we should tell that the cascade at a fc-step, or the trajectory 
at a fc-point, is locally stable with respect to the initial point ip. Similarly, by means of 
the ultralocal multiplier 

My?) = 7 — = r~\ ( 20 ) 

one may characterize the local stability with respect to the variation of the value yu-i 
preceeding the point yk, which takes place for 

\M<P)\<1- (21) 

For brevity, we shall say that the cascade at a A;-step is locally stable if inequality (19) 
is satisfied and it is ultralocally stable when condition (21) is valid. The latter can also 
be called the contraction condition. The images of the corresponding multipliers for the 
domain of the variable x are defined as 

m k (x) = /j, k (F (x,u k (x))) , m k (x) = JI k (F (x,u k (x))) , (22) 
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in accordance with Eq. (4). Note that the local multipliers (18) define the local Lyapunov 
exponents 

X k (<p) = ^\n\pi k (<p)\ (A; = 1,2,...), (23) 

so that the stability condition (19) is equivalent to X k ((p) < 0. Some more information on 
dynamical systems and their stability can be found in Refs. [47-50]. 

The local stability of an approximation cascade {y k }, that is, the local stability of its 
trajectory {y k (<p)}, means the local convergence of the approximation sequence {/&(#)}, 
bijective to this trajectory. In other words, if Eq. (19) is correct, then the /c-order approx- 
imation is closer to a fixed point, that is more accurate, than the initial approximation. 
And when Eq. (21) is valid, the fc-order approximation is better than the (k — l)-order 
approximation. The stability conditions (19) and (21) are sufficient for local convergence, 
but not necessary. So, when such conditions do not hold, this does not necessarily mean 
local divergence, but only tells that convergence cannot be guaranteed. 

In order to characterize the type of convergence, let us consider the variation 

tiyM = yk(p + Sip) - y k (<p) = Hk(<p)5<p . 

Because of the relation (23) between the local multipliers and Lyapunov exponents, we 
have 

\8y k (ip)\ = \8ip\ exp{X k {p)k} . (24) 

For the stable process, when \ k (<p) < 0, the value \5y k (ip)\ exponentially decreases with 
k. Therefore one tells that the calculational procedure is exponentially stable [51]. This 
is equivalent to saying that the approximation sequence {fk(x)} exponentially converges, 
since this sequence is, by construction, bijective to the approximation cascade trajectory 
{y k (f)}- Similarly, we could consider the variation 

5yl( V ) = yl( V + Sip) - y* k {p) = vL k (tp)6<p , 

which results in the expression 



|^)|<C|M exp{\ k (p)k}, 
12 



(25) 



provided that \5yl((p)/6yk(<p)\ < C. The latter inequality holds true for the stable 
process, when \k(f) < 0, then Vk(ip) — > 0, as k — > oo, and ~~ Uki^)-* so that 



When the values of the sought function f(x) can be found numerically, the accuracy 
of our calculational procedure can be characterized by the relative errors 



Analysing these, we can judge whether the accuracy of the procedure becomes higher with 
increasing approximation order k. But if f(x) is not known, how could we decide about 
the validity of approximations? In such a case, we have to be guided by the stability 
analysis dictating the following strategy: (i) When both stability conditions (19) and (21) 
hold true, then we can be certain that we are approaching the correct answer. In this case, 
the minimal time r in the evolution integral (17) corresponds to the minimal number of 
steps, that is to r = 1, necessary for reaching starting from (ii) When condition 
(19) holds but (21) does not, this means that the A;-order approximation is certainly 
better than the initial approximation but may be worse than the preceeding (k — 1)- 
order approximation. In this situation, the accuracy of the sought approximation can be 
improved by employing the middle-point method, well known in numerical analysis [51]. 
The middle point can be defined in two ways, explicit and implicit. For the evolution 
integral (17), the explicit definition of the middle point implies that we have to associate 
fl +1 with the middle of the interval between k and k + 1, that is, we have to set r = \. 
The implicit definition of the middle point between two approximations, say -Ffe+i and 
Fk, is to take the average ^(Fk+i + Fk) of these approximations [51]. Using this average, 
instead of F k+ i, in the cascade velocity (16) yields the factor | for the latter. This 
evidently is identical to putting r = \ in the evolution integral (17). So, the explicit and 
implicit ways of defining the middle point are equivalent to each other, (iii) Finally, if 
both conditions (19) and (21) are not valid, we cannot trust to the result of calculations. 
This can happen if the choice of control functions or the initial approximation, or both of 
them are not appropriate. Then we need to change these so that to achieve the stability 



| | -> 1. 
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of the procedure. 

One more question could be as follows. Assume that we can realize our calculational 
procedure in two ways, say using different initial approximations. Then which of these 
ways should we prefer? The answer to this question straightforwardly follows from the 
stability analysis: That of several possible ways is preferable which is the most stable. 

Thus, in the method of self-similar approximants, there is a unique opportunity to 
control the validity of calculations not knowing the sought function and having a small 
number of terms. In order to unambiguously prove all this, we consider in the following 
sections several examples for which numerical results are available. Comparing the latter 
with our predictions, everyone can directly observe that the method does work. 

3 Effective Potentials 

Consider the case of the so-called zero-dimensional <p 4 theory. The effective potential, or 
generating functional, has the structure 

/(<?) = - In Z(g) , Z(g) = ±= J exp {-S(g, <p)} dip , (26) 

with the action 

S(g,tp) = tp 2 + g<p\ (27) 

where the field ip is a real variable, f G R, and the coupling parameter g e [0, oo). 
Expanding f(g) in powers of g yields, as is well known, divergent series. 
To introduce control functions, we rewrite the action (27) as 

S(g, <p) = S (<p, u) + AS(g, <p, u) , 

S (<p, u) = M V , &S(g, <p, u) = (1 - « V + 9^ ■ (28) 
When AS* = 0, we have, in the place of Eq. (26), 

F (g,u) = lnu. (29) 
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Invoking the expansion in powers of AS, we find 



with the notation 



and the coefficients 



k 

F k (g,u) = F^fau) + £ c kp a k ~ p (3 p , (30) 

p=0 



I 3o 
« = = l-_, /3 = /3(u) = ^ (31) 



c io — — 2 ' Cn — 4 ' ° 20 ~~ — 4 ' C21 ~~ 2 ' ° 22 ~~ _ 3 ' 
I _ 3 4 II 

C30 = — g , c 31 — ^ , C 32 — — - , C33 — — , 

I 10 11 34 

Qo — — g , c 4 i — 1 , c 42 — — — , c 43 — — , c 44 — — — , 

and so on. As is seen, expression (30) is not a series in powers of g but a bipolynomial 
with respect to the variables a and (5 defined in Eq. (31). 

Control functions Uk{g) are given by the fixed-point condition (15). This means that 
for odd k — 2n + 1, we use the equation dFk/duk = 0, while for even k = 2n, for which 
the latter equation has no solution, we put u k = Uh-i- This procedure gives 



Uk{g) 



±(l + y/l + 12s k g 



1/2 



(32) 



where Si = S2 = 1, S3 = S4 = 2.239674 , . . .. For large k ^> 1, one has s k ~ k £ , with 
1 < e < 2. Substituting u k (g) into Eq. (31), we have the relation 

®(u k (g)) = a k (g) = s k (3(u k {g)) , (33) 

in which 

i-yi + i2s fe fl 

«fc 5) = 1 + 7 • 34 

6s k g 

To construct the approximation cascade, we directly follow the procedure of Sec. 2. 
The sole unimportant difference is that, instead of the variable x there, we have here the 
variable g. The equation (4) for the expansion function now reads lnu k (g k ((p)) = if, which 
yields 

9k(<p) = 3^ {e 2ip ~ l) • (35) 
15 



The transformation (5) takes the form 

k 

y k (ip) = ^ + ^A kp a p {^>) , 
P =i 

where the coefficients are 

4*=E^ (p<*) 

m=0 *fe 

and the notation 

«(</?) = a k (g k ((p)) = 1 - e~ 2ip 
is used. All coefficients (37) are negative, e.g. 

Mx = ~\ , An = A u , A 22 = , 

A 31 = -0.388377 , A 32 = -0.093205 , A 33 = -0.016011 , 
A 41 = A 31 , A 42 = A 32 , A 43 = A 33 , A 44 = -0.003606 . 
The cascade velocity (16) is 

v k ((p) = B k a k+1 ((f) , B k = A k+1M+1 . 

Invoking the notation 

f* k (g) = In y/l + z* k (g) , f k (g) ee In Jl + z k (g) , 
we can write the evolution integral (17) as 

r*i + i (i + z) k J 

/ zk+l dz = 2B k r , 

where z k = z k (g) and z k = z k {g). Integrating Eq. (41), we obtain the equation 

4+1 = 2t exp{p(-i-)-p(i) +2Bi ,}, 
defining z% +1 (g), where 

k-l % k-p ^| 
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To check the stability of the approximation cascade {yk}, we calculate the local mul- 
tipliers (18). This gives 

k 

M = 1 + 2[1 - a( V )]J2 P A kp a p ~\ V ) . (43) 

P =i 

It is interesting to note the relation between the subsequent multipliers, 

Pk+i(<p) = M^) + A k(<f) , (44) 

where 

Ak(ip) = 2(k + 1) B k [1 - a(<p)] a k {ip) . (45) 

From formulas (35) and (38), it follows that the domain of g G [0, oo) corresponds to 
tp G [0, oo) and a((p) G [0, 1]; that is, the effective expansion parameter a(ip) in Eqs. (36) 
and (43) is always less than 1 for arbitrary g G [0, oo). Also, for all ip G [0, oo) the stability 
conditions (19) and (21) are valid. Therefore, we put r = 1 in Eq. (42). 

It is convenient to characterize the accuracy of the approximants (40) by the corre- 
sponding maximal errors e* k = sup s |££(<7)|, Sk = sup 9 \sk(g)\- Comparing these approx- 
imants with direct numerical calculations for the function (26), we have S\ = 7%, e 2 = 
4%, e 3 = 0.2%, e 4 = 0.2% . . ., and e* 2 = 3%, e* 3 = 2%, el = 0.1% . . ., which demonstrates 
good local convergence of the sequence {fk(g)} of the self-similar approximants . 



4 Anharmonic Oscillators 

In + 1 dimensions the ip A theory reduces to the familiar anharmonic oscillator prob- 
lem. Therefore, it is instructive to illustrate the approach for the eigenvalue problem of 
quantum mechanics. Let us consider the Hamiltonian 

ffW = -2^ + ^' < 46 » 
in which m,A,v > 0, and x G R. This Hamiltonian, by scaling its variable, can be 
presented in the dimensionless form 

H = -\^ +g *. (47) 
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To return from Eq. (47) to Eq. (46), one has to make the substitution 

H(t\ A 2 

H -> tl±±i x ^ Xj g^l, uj 2 +» EE — . (48) 

So, we need to find the eigenvalues e(n, g), where n — 0, 1, 2, . . ., of the Hamiltonian (47), 
and then, according to Eq. (48), we have to put g — 1 obtaining the sought spectrum 
e(n) ee e(n, 1). The latter can be compared with the direct numerical calculation of the 
Schrodinger equation, as well as with the quasiclassical approximation [52] which gives 



7T / 1\ / iA „ /2 + iA ,„ (\ 



2v/ (2+i/) 



(49) 



where T(-) is the Gamma-function. 

For v — 2, the Hamiltonian (47) is that of harmonic oscillator. This suggests to take 
for the initial approximation the harmonic oscillator 

1 d 2 u 2 9 . . 

H ° = ~2d7' + T ' (50) 

in which u is a trial parameter that will generate later control functions Uk{g). The 
eigenvalues and eigenfunctions of H are, respectively, 

1\ _ L (m/tt) 1 / 4 

(2™n!) ] 

where n = 0,1,2,... and H n {-) is a Hermit polynomial. By means of the Rayleigh- 
Schrodinger perturbation theory, with respect to the perturbation H — H , one can find 
the sequence {E k (n, g, u)} of the /c-order approximations for the eigenvalues of the Hamil- 
tonian (47). For what follows, it is convenient to introduce the notation 

E k (n,g,u) ee (n + -] F k (n,g,u) (51) 



E (n, g,u)= \n+-j U) ip n (x) = , 2 „ 1/2 exp \-- x j H n [^/u x J , 



2 

and to consider the sequence {F k (n, g, u)}, with k — 0, 1, 2, . . . Defining control functions 
u k (n,g) from the fixed-point equation (15) and substituting them into Eq. (51), we get 

e k (n) = (n + ^ lim F k (n, g, u k (n, g)) . (52) 

And from the evolution integral (17), we find the self-similar approximant e* k {n). 
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Quartic oscillator. Let us realize the above program for v = 4 in the Hamiltonian 
(46). For the perturbative terms F k) defined in Eq. (51), limiting ourselves by the second- 
order approximation, we have 

F (n, g,u) = u, F^n, g, u) = F (n, g, u) - - + ^ j n , 

F 2 {n,g,u) = Fx{n,g,u) - | + ^ 7n - ^ a n , (53) 

where 

n 2 + n + ± 1 1 27 / 1\ / 1 

For the control function we obtain Ui(n, g) = (67„ sO 1 / 3 . Equation (4) now reads 
Uk(n, gk(n,ip)) = if, which results in the coupling function gi(n,<p) = y? 3 /67 n . The endo- 
morphism (5) is written as yk(n,ip) = /ik(n)(p, where 

Mn) = l, ^(n) = \ + A(n) , A(n) = ± - ^ . (54) 

The cascade velocity (16) becomes v\{n,<p) = A(n)(p. Formula (52) gives 

3 / 1 



ei (n) = - \n + - j (6 7n ) 1/3 , e 2 (n) = 



1 + |A(n) 



ei (n) . (55) 



And from the evolution integral (17), we obtain 

e^{n) — ei(n) exp {A(n)r} . (56) 

To decide what effective time r to choose in the self-similar approximant (56), we 
have to analize the multipliers (18) and (20). For the local multipliers (18), we get 
/ik(n, V 9 ) = / i fe( n )j with /ifc(n) given in Eq. (54). When varying n — 0, 1, 2, ... 00, we have 

35 . . 109 1 A/ . 1 

— < Pain) < , < A(n) < . 

4g _ ) - M4 , 4g - \ ) — 144 

Defining the maximal multipliers 

^ fc = sup \nk{n)\ , ~p k = sup|7Z fc (n)| , (57) 
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n n 



we obtain \i\ = 0.75, /i 2 = 0.756944, /Z 2 = 1.009259. It is the maximal multipliers (57) 
that are to be considered in the stability analysis. This is because the sequence {F k }, 
as in Eq. (53), has been derived employing the Raleigh-Schrodinger perturbation theory, 
which involves the summation over quantum numbers. Therefore, the terms F k (n, g,u), 
and consequently the energies (51) and (52), with different quantum numbers n, cannot 
be treated as completely independent. 

Since ii k < 1, the procedure is locally stable; but ~p 2 is slightly larger than 1, so that 
there is a week ultralocal instability. Thus, we need to take the middle point r = |; 
although, because the ultralocal instability is so weak, we could expect that the results 
should not essentially differ between the cases r = \ and r = 1. 

To check the prediction of the stability analysis, we may compare the energy levels (55) 
and (56) with the results of numerical solution of the corresponding Schrodinger equation 
[53]. Varying the quantum numbers n — 0, 1, 2, . . ., we define the maximal errors 

£wkb = sup \e WKB (n)\ , e k = sup \e k (n) | , e* k = sup \e* k (n)\ . (58) 

n n n 

For the latter we find e WKB = 1.5%, E\ = 2%, e 2 = 0.8%, e 2 = 0.4%. The error e* 2 does 
not change for the alternative choices of r = 1 and r —\. 

Sextic oscillator. For the power v = 6 in the Hamiltonian (46), we have the sequence 
F (n, g,u)=u , F^n, g, u) = F (n, g, u) - - + — n n , 

F 2 (n, g, u) = Fi(n, g,u) - - + K n - ^7 A» , (59) 
which is defined in Eq. (51), and where 

Kn = n 2 + n + - , (3 n = 786n 4 + 1572n 3 + 5324n 2 + 4538n + 3495 . 

Then we follow the same steps as for the quartic oscillator. We find the control function 
Ui(n,g) = (15«; n g) 1 ^ and the coupling function gi(n,ip) = (p 4 /lhn n . Constructing the 
approximation cascade, we come to the same endomorphism y k (n,ip) = /i k (n)ip but with 
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From Eq. (52), we get 

2 / I 



2 / 1 \ 

e i( n ) = 3 [n + - J (15/t n ) 1/4 , e 2 (ra) = 



1 + ^A(n) 



ei(n) . (61) 



The integral (17) gives the same expression (56), but with A(ra) defined in Eq. (60). 

For the local multipliers (18), we have the same form as in the previous case, with 

Hk( n ) given in Eq. (60). Varying the quantum numbers n — 0, 1, 2, . . ., we find 

311 , . 273 19 A , . 49 

< Mn) < 77^ > ~T7^ < A ( n ) < 



540 ~ ^ y 1 ~ 400 ' 1200 ~ v ' ~ 540 

The maximal multipliers (57) are /ii = 0.6667, /i 2 = 0.6825, JJ 2 = 1.023699. Again we 
see that the procedure is locally stable since fij, < 1, but with a weak ultralocal instability 
because of ~p 2 > 1. Therefore, we again have to take in Eq. (58) the middle point 
t — \. For the maximal errors (58), with respect to numerical calculations [54], we find 
£wkb = 4%, S\ = 7%, e 2 = 8%, e* 2 = 2%. Taking r = 1 slightly increase the error to 
3%. Note that in all cases the maximal error occurs at the ground-state level. 

Octic oscillator. With the power v = 8 in the Hamiltonian (46), we find the sequence 
of perturbative terms 

F (n,g,u) = u , Fi(n,0,u) = F (n,g,u) - - + — a n , 

F 2 (n, «) = Fi(n, g, u) - - + — a n - S n , (62) 

where 

_ n 4 + 2n 3 + 5n 2 + An + 3/2 
n + 1/2 

<5 n = 3985n 6 + 11955n 5 + 74904n 4 + 129883n 3 + 277901n 2 + 214952n + 135030 . 

The control function is Ui(n, g) = (35<x n g) 1 ^ 5 and the coupling function becomes gi(n, ip) = 
(p 5 /35<j n . The approximation cascade is described by the same endomorphism with 

*(») = §■ * = f + A<»), Al,,),^. (63) 



n 



For varying n — 0, 1, 2, . . ., one has 



4319 , N 5083 3031 . . . 183 

< »2(n) < — - , -tt^tt < A(n) < 



11760 ~ ^ v 7 ~ 7840 ' 11760 ~ v ' ~ 7840 
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The energies (52) are 

5 / 1 



e iW = g ( n + 2) ( 35(T ™) 1/5 ' e2 ( n ) = 



1 + -A(n) 

5 



ei(n) . (64) 



And the self-similar approximant e^n) again has the same form (56). 

The local multipliers are given by /x/c(n). For the maximal multipliers (57) we find 
lii = 0.625, /i 2 = 0.648342, /Z 2 = 1.037347. This shows that the procedure is locally 
stable but is not ultralocally stable. So, the self-similar approximant (56) has to be 
defined at the middle point r = \. 

For the maximal errors (58), with respect to numerical calculations [54], we obtain 
£\vkb = 7%, E\ = 13%, e 2 = 34%, = 3%. The accuracy of the self-similar approx- 
imant (56) for t = I is essentially higher than that for r = 1, for which case the error 
reaches 13%. 

Multidimensional oscillator. In the previous three cases we have shown that the 
method works well for different one-dimensional anharmonic oscillators. Now we want to 
demonstrate that it is equally applicable for anharmonic oscillators of arbitrary dimen- 
sionality. For this purpose, we consider the spherically symmetric quartic oscillator in 
.D-dimensional real space. The corresponding radial Hamiltonian reads 

where to, A > 0; I — 0, 1, 2, . . . ; r > 0, and D — 1, 2, 3, ... is a real-space dimensionality. 
This Hamiltonian, by the appropriate scaling, can be reduced to the form 
„ Id 2 (2/ + D-3)(2/ + D-l) 4 

H = - 2 d^ + ^ +9r ■ (66) 

One can return from Eq. (66) back to Eq. (65) by means of the substitution 

rr H{r) MV /3 

H — > , r — > vto r , g — > 1 , lo _ 



uj \m? 
We start perturbation theory with the harmonic D-dimensional oscillator defined by 
the Hamiltonian 

1 d 2 (2l + D-3)(2l + D-l) u 2 2 
H °-~2d^ + 8^ + T r ' 
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in which u is a trial parameter. The corresponding eigenvalues are 



E (n, Z, u) = (^2n + I + — j u , 



where n — 0, 1, 2, . . ., and the eigenfunctions are 

2n \ u i+D/2 nV2 



Xm(r) 



r l2l + D-l)/2 exp ( _U r2 j L ( 2i+ D-2)/2 (w , 2) ? 



_r(n + / + D/2)_ 

where ^(r) is the associate Laguerre polynomial. The latter can be presented in several 
forms, so we write down below that one we use in what follows: 

j) (r) = L e r r -« *L ( e -r r n +l) = ^ T(n + Z + l)(-r) m 

" v ; n! dr n v ' ^ T(m + Z + l)m!(n - to)! ' 

For the terms of perturbation theory, we use the notation 

E k {n, Z, g, u)=(2n + l + F k (n, Z, g, u) . (67) 

Then, in the second order we have 

F (n,l,g,u) = u , Fi(n,Z,#,u) = F (n,l,g,u) - | + ^ -y nl , 

F 2 (n, Z, g, u) = F^n, Z, g, u) - | + ^ 7„/ - ^ a nZ , (68) 

where 

_ D (l + D/2)(l-2 + D/2) 

^ = 2n + l + -- 3{2n + l + D/2) , 

27/ , D\ ( , L>\ 2 

a nl = 1 + — (^2n + / + —J lnl -{2n + l + -) . 

The one-dimensional case can be recovered with the change D — > 1, 2n — > n, Z — > 0. 

All steps of the procedure are the same as earlier, so we do not go into much details. 
For the control functions we find Ui(n, Z, g) = (G^fnig) 1 ^ 3 , and for the coupling function we 
get 0i(n,Z,<p) = y? 3 / 6 7n2- 

The approximation cascade is characterized by the endomorphism yk(n, Z, ip) = //^(n, Z)</?, 
in which 

/i!(n,Z) = ^, /i 2 (n,Z) = ^ + A(n,Z) , A(ra, Z) = ^ - , 
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and the cascade velocity is vi(n, /, (p) = A(n, l)ip. 
Similarly to Eq. (52), we define 

e k (n, l) = ^2n + l + lim F k (n, I, g, u k (n, I, g)) . 

This yields 



l + ^A(n,/) 



e 1 (n,l). (69) 



ei (n, I) = - (2n + I + -) (6 7ni ) 1/3 , e 2 (ra, /) = 

And the evolution integral (17) gives us the self-similar approximant 

e*(n, = ei(n, /) exp {A(n, Z)r} . (70) 

The stability analysis here is practically the same as for the one-dimensional quartic 
oscillator, with the same conclusion that the procedure is locally stable but with a very 
weak ultralocal instability. Consequently, in the approximant (70) one has to take the 
middle point r = |, although the results should not be much different from the case r = 1. 
The accuracy should increase with the increasing dimensionality D, since increasing D 
supresses ultralocal instability. For example, when D — > oo, then 7^ — > D/3, a n \ — > D 2 /5, 
and A(n,l) — > 0, so that /Z 2 — > 1. These conclusions are in perfect agreement with 
calculations. Thus, the maximal error of the self-similar approximant (70) is e\ = 0.4% 
for D = 1 and it is ^ = 0.3% for D = 3. 



5 Arbitrary Powers 

In the previous section, we considered several Hamiltonians with potentials having even 
integer powers. However, this restriction is not principal, and the method can be applied to 
the case of potentials with arbitrary, noninteger as well as integer, powers. To demonstrate 
this, let us consider a spherically symmetric three-dimensional problem with the radial 
Hamiltonian 
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in which m, A > 0; / = 0, 1, 2, . . . ; r > 0; and v > is an arbitrary positive number. By 
scaling, it is again convenient to reduce this Hamiltonian to the simplified form 



The return from Eq. (72) to (71) is made by means of the substitution 

H^^-, r^V^r, <?-l, cu 2 ^ = — . 

Hamiltonians of this type have been considered in a number of papers (see e.g. [31,55-58] 
and references therein). The importance of potentials with various powers is due to the 
fact that such potentials model well confining forces between heavy quarks and lead to a 
reasonable description of the quarkonium spectrum [59]. The quasiclassical approximation 
gives [59] the eigenvalues of the Hamiltonian (72) in the form 



ewKB(n,l) = 



7f / , 3\ / v\ „ /2 + i/\ ,„ /l 



- 2n + i + - i + _ r — /r - 



2j//(2+i/) 



(73) 



Applying our method to the eigenvalue problem for the Hamiltonian (72), we start 
with the harmonic Hamiltonian 



1 d 2 1(1 + 1) u* 



whose eigenvalues are 

E (n,l,u) = (2n + l + ^ju . 
Using the Rayleigh-Schrodinger perturbation theory, we can find the /c-order terms 



E k (n, I, g, u) = ( 2n + I + - ) F k (n, I, g, u) . (74) 



2 

Then, defining control functions Uk(n,l,g), we come to the renormalized eigenvalues 

e k (n, l)= (2n + l + lim F k (n, I, g, u k (n, I, g)) . (75) 
In the second order, we obtain the sequence 

F (n, l,g,u) = u, F^n, /, g, u) = F (n, l,g,u)-^ + ^ A nl (^j , 

25 



F 2 (n,l,g,u)^F 1 (n,l,g,u)-^ + ^^B nl (^j -JL-Crafy , (76) 



with the coefficients 



A n ,(u) = 



(2n + / + 3/2)r(n + / + 3/2) 
C nl {y) 



(2n + / + 3/2)r(n + / + 3/2) ' 

r '« + i)t + iM-(" + ' + ^)Ci(") 



E 







2 


p! 


4»" 





(2n + Z + 3/2)r(n + Z + 3/2) p=Q ^ n) (p - n) F(p + Z + 3/2) ' 
Here we use the notation 



JO 



(77) 



where 14(0 is an associate Laguerre polynomial. The way of calculating this integral is 
explained in Appendix A. As is seen from the sequence (76), we have there an expression 
in powers of g/u l+v l 2 . 

Following the same steps as in the previous section, we find the control and coupling 
functions 



ui(n,l,g) 



vA nl \-)g 



2/(2+1/) 



9i{n,l,<f) 



1+I//2 



vA nl {vj2) ' 

The approximation cascade is defined by the endomorphism yk(n, Z, ip) = Hk( n , Ov 9 with 

2 + v 



//j (//./) = z -^- , i^ 2 (n, I) = + A(n, Z) , 

5 ni (z//2) 0/2) 1 



A(n,Z) = 



2uA nl (u/2) 2v 2 A nl (v/2) 8 



(78) 



For the values (75), we have 



ei(n,Z) = (2n + Z + -^ 



2/ 2i/ 



2/(2+1/) 



2 is . . 7N 



2, 

ei(n,Z) . 



e 2 (n, Z) = 

The evolution integral (17) yields the self-similar approximant 

el(n, Z) = ei(n, Z) exp {A(n, Z)r} . 



(79) 



(80) 
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The local multipliers (18) are I, ip) = /ifc(n, /), with the right-hand side defined 
in Eq. (78). The analysis of these multipliers shows that the stability properties depend 
on the power v of the potential in the Hamiltonian (72). For v ^ 2, the procedure is 
either locally unstable but then it is ultralocally stable, or the procedure is locally stable 
but ultralocally unstable. In both these cases, we have to take in the approximant (80) 
the middle point r = \. The case v = 2 is special. Then all multipliers ji k = JI k = \, 
which corresponds to the neutral stability. This means that we are already at a fixed 
point of the cascade trajectory, and there is no motion any more. Really, our calculations 
show that for v = 2 and any r, we have for all eigenvalues ei(n, /) = e 2 (n, /) = e^n, /), 
which coincides with the exact solution for the harmonic oscillator. To illustrate what 
are the values of typical errors characterizing the accuracy of different approximations, 
we present in Table I the results of calculations of the ground-state energy for various 
powers v. We give there the errors E\ and e 2 for the renormalized expressions from 
Eq. (79), the error e* 2 for the self-similar approximant (80), and also the error of the 
quasiclassical approximation (73). These errors are computed by comparing the results of 
our calculations with accurate numerical data [55,57,59] obtained by the direct numerical 
solution of the Schrodinger equation. The results for the approximant (80) correspond 
to the middle-point case r = |, which is taken according to the stability analysis. For 
comparison, we also made calculations setting r = 1. The latter case leads to the errors 
that are close to e\ from Table I, when v < 6. However, for v > 6, the errors quickly 
increase becoming, e.g., —6% for v = 8 and —19% for v = 10. Hence, the middle point 
choice r = I is especially preferable for large powers of the potential. 



6 Initial Approximations 

Now we aim at analysing the following problem. Assume that we can take not just 
one but two or more initial approximations that can be employed as starting points for 
constructing the corresponding approximation cascades. Then how to decide which of 
these initial approximations is preferable? We show that the stability analysis again can 
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give us an answer to this question. 

For the purpose under consideration, it is convenient to take the logarithmic potential 
whose behaviour is such that it can be modelled, to some extent, either by the harmonic 
oscillator or by the Coulomb potential. We mean here the three-dimensional spherically 
symmetric problem with the radial Hamiltonian 

ff M = -^ + ^ + Bh i- (81) 

in which m,B,b > 0, and r > 0. As usual, we transform the Hamiltonian (81) to the 
scaled form 

1 d 2 1(1 + 1) 

The return from Eq. (82) to Eq. (81) is made by means of the substitution 

H(r) r B 1 

H -> , r -> - , 9^ — , cj- 



u b to mb 2 

Note that the logarithmic potential is also often used for describing confining forces be- 
tween quarks [59] , and the quasiclassical approximation for the spectrum of the Hamilto- 
nian (82) is known to be 



ewKB(n,l,g) = |ln 



7T / ,3 



, 2n + / + 
2g V 2 



(83) 



Harmonic potential. Let us choose for the initial approximation the harmonic- 
oscillator Hamiltonian. Then for the function 

j-i/ r x_ E k {n,l,g,u) 

F k (n,l,q,u) = -, -— 84 

v ' ' y ' ! (2n + / + 3/2) v ; 

we derive the sequence of perturbative terms 

F Q (n,l,g,u) = u , 

. u ip(n + / + 3/2) — In if 
(n ' Z ' 9 > U) = Fo(n ' l >9> U) - 2 + 9 2(2n + / + 3/2) ' (85) 

9 

F 2 (n, I, g, u) = Fl (n, l,g,u)-- + ^ + f + - - C nl , 
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where 



a 



nl 



^( x ) = ^ In r(x) , 



E 



(p V 



(2n + l + 3/2)r(n + I + 3/2) p=Q ^ n) (p - n)r(p + Z + 3/2) ' 
Here 7^ s denotes the integral 



4- f> 1/2 e" 4 lntL^/ 2(t) 4 + i/2 (0 

JO 



(86) 



The properties of this integral are described in Appendix B. 

Following the standard scheme of our method, we find the control and coupling func- 
tions 



ui(n,l,g) = 



9 



9i(n,l,ip) = (2n + l + up . 



2n + I + 3/2 ' 

Then we construct the approximation cascade {y k }, with the trajectory points 



yi(n,l,(p) = 



I - \n(p + ip [n + l + 



y 2 (n,l,(p)=y 1 (n,l,<p)+vi(n,l,<p) , 

(87) 

and with the cascade velocity vi(n, I, (p) = A(n, l)ip, where 

A(n, I) = 1 - 1 (2n + I + C nl , A(0, 0) = -0.018288 . 
For the renormalized expressions 

e k (n, I, g)= (2n + I + F k (n, I, g, u k {n, I, g)) , 



(88) 



we have 



ei(n,l,g) 



9 



1 + In '— +ip[n + l + 

9 V 



e 2 (n, I, g) = e^n, I, g) + A(n, Z) 9 . 
From the evolution integral (17), we obtain the self-similar approximant 



(89) 



el(n,l,g) = e x {n, Z, g) exp {A(n, Z)r} 



(90) 
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The qualitative behaviour of the energy (90) is as follows. At small g — > +0, this energy- 
is positive and tends to zero as —glng. As g increases, e\ increases till the maximum 
value 

e* 2 (n,l,g max ) = ^ g max (n, I) exp{A(n, l)r} , 

occurring at 

9ma X (n, l) = [2n + l + -) exp |^ (n + I + -) J . (91) 

With the further increase of g, function (90) decreases to zero, e^n, l,g ) = 0, at the 
value 

g (n,l) = eg max (n,l) , (92) 
after which e\ remains negative. Using notation (92), we may present Eq. (90) as 

e*(n, l,g) = 9 - In (^j exp {A(n, l)r) , (93) 

where g = g (n,l). 

To estimate the values of the characteristic couplings (91) and (92), consider the case 
n=l=0, when 

9maM 0) = l e 2 " c = 1.555746 , g (0, 0) = jj e 3 ~ c = 4.228956 . (94) 

where C is the Euler constant. 

In the opposite case, when either n or / is large, since ip(x) ~ lnx, as x — > oo, we have 



gmax I) — * 



(95) 



2n 2 (n — > oo, I < oo) 
I 2 (n < oo, I — > oo) . 

The local multipliers (18), for the trajectories points (87), are 

l*i(n,l,<p) = 7^ ( n + l + ^) - g ' ^(n,l,ip) = /Xi(n,Z,</?) + A(n, Z) . (96) 

The images of the multipliers (96), defined in Eq. (22), become 

, , , 1 . / , 3\ 1 . 2n + / + 3/2 
mi(n, /, = - ^ I n + / + - J + - in 



2 r V 2/2 5 

m 2 (n,l,g)=m 1 (n,l,g) + A(n,l). (97) 
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With definition (91), we may write 

m 1 (n,l,g) = 1 -\n 9max{nJ) . (98) 
^ 9 

Varying the quantum numbers n and /, one can check that A(n, /) is always negative and 
the ultralocal multiplier m 2 (n, l,g) = 1112(11, 1, g)/mi(n,l, g), which is the image of the 
multiplier (20), is less than one for all 11, I, and g. That is, the procedure is ultralocally 
stable. However, it is locally stable not for all coupling parameters but only for those 
satisfying the inequalities 

- 2 < < e 2 , (99) 

6 g-max 

where g max = g m ax(n, I) is given by Eq. (91). For instance, if n — I — 0, then the interval 
(99) is 

0.210547 < g < 11.495494 . (100) 

In the region of stability (100), we may put r = 1. For illustration, we present in Table 
II the results of calculations for n = and I — 0, 1, 2, 3, 4, and for the coupling parameter 
g = 0.5. We show the percentage errors E\ and e 2 for the renormalized expressions (89), 
the error e* 2 of the self-similar approximant (90), and the error Swkb of the quasiclassical 
approximation (83). The comparison is made with respect to numerical data [55,59]. 

Coulomb potential. Now let us turn to the case when the Hamiltonian 

h - 1 d2 1 ^ + 1 u non 

Ho --2d^ + ~2^ + r (101) 
with the Coulomb potential is chosen for the initial approximation. The corresponding 
eigenvalues and eigenfunctions are 

u 2 

E (n,l,g,u) = 



2(n + l + l) 2 ' 



Xm(r) 



ni u 



(n + 2l + l)\ 



1 ' 1 ' 2ur \ 1+1 ( ur \ , 2J+1 / 2n> 



ex P -Z-TTT K 



n + l + \\n + l + \) V n + l + l) n \n + 1 + 1 



In calculating matrix elements, we meet the integral 

JL = [°°t 2l+2 e" 4 \nt L^(t) Lf+\t) dt , (102) 
J 
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whose properties are specified in Appendix C. The whole calculational procedure is very 
similar to the previous subsection, because of which we shorten here technical details. 
This is especially justified by the fact that the procedure starting with the initial approx- 
imation corresponding to the Hamiltonian (101), as will be shown, is less stable than that 
beginning with the harmonic-oscillator Hamiltonian. 
The spectrum in the first approximation reads 

u 2 ( n + 1 + 1 \ 

E ' ("'''»' "> = 2 ( n + ; + 1)2 + g ¥a + D "') ■ (103) 



where 



2n + 1 

p "" 2faT7TT) + ^ ( " + 2 ' + 2) - 



For the control and coupling functions, we find 

ui(n,l,g) = (n + I + l)y/g , g x {tp) = -2tp . 

Then the point of the approximation cascade, bijective to the approximation (103), is 

yi (n, I, ip) = -(p[l - ln(-8</?) + 2D nl ] . 

Substituting the control function into spectrum (103) gives 

ei (n,l,g) = ^(l + \nj- + 2D n )j . (104) 

The self-similar approximant e* 2 acquires the form of Eq. (90). 

The overall qualitative behaviour of the spectrum is the same as in the previous 
subsection. For g — * 0, the energy tends to zero as —glng. With increasing g, the energy 
reaches a maximum at the point 

g ma x(n, I) = J exp(2D ni ) . (105) 

Then the energy diminishes to zero at the value g = g (n,l) having the same relation 
go — zg m ax, as in Eq. (92). For the ground-state level, we now have 

g max {0, 0) = - e 3 - 2C = 1.582925 , g (0, 0) = - e 4 " 2C = 4.302836 . (106) 
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The energy (104) can be written as 

ei(n, «,(/) = | ln(j^ , (107) 

where g = go(n, I), and e\ can be presented as in Eq. (93). 

The characteristic values of g max , given by Eqs. (105) and (91), are very close to each 
other. This is seen from Eqs. (94) and (106), as well as comparing Eq. (95) with the 
corresponding behaviour of Eq. (105) yielding 



(e 2 /A)n 2 (n — >■ oo, I < oo) 
I 2 (n < oo, I — > oo) . 



9max iP'i — * 

For the local multiplier (18) and its image (22), we find 



Pi(n,l, <p) = ln(-8</?) - 2D nl , mi(n, Z, #) = In — -. (108) 

9max [f^i «) 

The stability condition \/ii(n,l,ip)\ < 1, \mi(n,l,ip)\ < 1 is satisfied in the region 

- < -9— < e , (109) 

6 9max 

where g max = g ma x(n,l) is defined in Eq. (105). For the ground state, the inequalities 
(109) give 

0.582326 < g < 4.302836 . (110) 

Comparing the stability properties of the calculational procedure for the considered 
cases, when either the harmonic oscillator or the Coulomb potential are used for the initial 
approximation, we come to the following conclusion. The region of stability, with respect 
to the coupling parameter g, is three times narrower in the latter case than in the former. 
Also, the values of the local multipliers in the latter case are about twice as large as for 
the former case. Consequently, the calculational procedure starting with the harmonic 
oscillator is more stable than that beginning with the Coulomb potential. This conclusion 
is in agreement with calculations showing that the first, more stable, case is also more 
accurate than the second, less stable; the error of the latter case being almost twice larger. 
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7 Conclusion 



The self-similar perturbation theory makes it possible to obtain quite accurate approxima- 
tions for those complicated physical problems for which only a few terms of a perturbative 
algorithm are available. The term "perturbation theory" here has to be understood in the 
general mathematical sense as a regular procedure yielding a set of subsequent approxi- 
mations. This approach does not require any small parameters of the physical problem 
considered. For instance, the coupling parameters can be arbitrary strong. Convergence 
is achieved by introducing control functions, by means of which the sequence of approxi- 
mations is reorganized to become convergent for any given parameters. 

The approach is based on the property of self-similarity between approximations. 
This property, for the family of endomorphisms representing the given approximation 
sequence, is a necessary condition for the motion in the vicinity of a fixed point. For the 
approximation sequence itself, this means a necessary condition for the fastest convergence 
[42]. 

Let us stress the difference between the self-similarity of approximations and the 
functional self-similarity which is the basis of any renormalization-group approach [46,60]. 
In the latter case, one looks for a symmetry of the considered function f(x) with respect 
to the scaling of the variable x, so that to find a relation f(Xx) = B(\, f(x)). From such 
a relation, it is easy to get the renormalization group equation xdf(x)/dx = f3(f(x)), 
where (3{f) = liniA->i dB(X, f )/d\ is called the renormalization group function. When 
the above relation is exact, one speaks of an exact renormalization group. However, in 
many cases such relations can be derived only approximately. As is clear from comparing 
this standard renormalization group technique with our approach, we do not consider 
the scaling of variables, but instead we are analysing a kind of scaling of approximation 
numbers. This what makes our approach principally different from any variant of the 
renormalization group techniques. 

The property of the self-similarity of approximations permits us to reformulate per- 
turbation theory to the language of dynamical theory and optimal control theory. Such 
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a reformulation allows us to resort to powerful methods of these theories in order to give 
a logical foundation for the whole approach. Thus, it immediately becomes clear that 
control functions are to be defined by the fixed-point equations. Then, considering the 
motion near a fixed point, we find nontrivial corrections essentially improving the ac- 
curacy of approximations. And, what is probably the most important, we can use the 
stability analysis in order to check the stability of the procedure and, respectively, the 
convergence of the approximation sequence. This is based on the analysis of the local 
multipliers that are directly related to the local Lyapunov exponents. 

What is especially valuable is that the stability analysis allows us to decide when we can 
trust to the results of calculations, even if we do not know exact answers. This also permits 
us do choose between several ways of calculations, differring by initial approximations. 
The general guide is always the same: We must select the most stable procedure. 

In the present paper, we have considered a variety of examples for which accurate nu- 
merical data are available. Then, directly comparing the latter with our results, we could 
explicitly demonstrate that the method does work. We think that such an explicit demon- 
stration is a necessary step before applying the method to more complicated problems for 
which no numerical data are known. 

Appendix A 

The integral (77) has the following properties that we have used in the process of 
calculations: 






71 ! 




where 



<Pn{x) 



T(n + x) 



x(l + x)(2 + x) . . . (n - 1 + x) , 



(p (x) = 1 , 



<pi(x) = x , 



x £ (— oo, +00) , 
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Also, 

ji m _L | j | 3\ r(n + / + 3/2) I> + j + 3/2) r(n + / + 5/2) 

Awt A J - ^ + « + g J ~\ °™ (n- 1)! s,n_1 n! ' n+1 ' 

In general, integral (77) can be presented as the sum 

, I> + / + 3/2) " (-iyrjp + l + 3/2 + u) 

nA ) s\ ^p!(n-p)!r(p + J + 3/2)^ P " J " 

Employing these properties, the coefficients in the sequence (76), for the case n — 0, can 
be simplified to 

. r(Z + 3/2 + 1/) , _ n 

Aw(Z/) " (Z + 3/2)r(Z + 3/2) ' " 1 ' 

" (Z + 3/2)r(Z + 3/2) ' B ° l{1) ~ 1 ' 



(Z + 3/2)r(Z + 3/2) ^ P p!r(p + / + 3/2) ' 
The value A(n, /) entering Eq. (78), for n — 0, writes as 



8 2z/ 2 ^ pp! T(p + Z + 3/2) 



Appendix B 



The integral (86), by employing the replica trick 

In t = hm 



can be expressed through the limit 

Jt = lim P nsM ~ 4,(0) 



V 



involving the integral (77). Taking this limit, we meet the expansion 

£^±^^l + ^> („-0), 
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in which 



d 00 ( 1 

= — lnT(x) = ^ 



q\P+1 p + x / 



A 1 



x 



\ p(p + x) X 



l -c, 



C = lim \ y - - Inn = 0.577216 



For n integer, one has 



V(n + 1) = £ - + ^(1), lKl) = -C, 
p=l p 



/ 1\ n 2 
*(n + -)=g— -2]n2-C, 



2 -In4-C = 0.036490 



In the course of taking the limit v — > 0, we also need to consider the corresponding 
asymptotic behaviour of the function y? n (— p — z/) defined in Appendix A. This asymptotic 
behaviour is 



<p n {-p-v) ~ < 



l)!p! 1/ (0<p<n-l) 
(-l) n n! {1 + [^(n + 1) - V(l)] v} (p = n), 

as 1/ — > 0. Using these equalities, for the integral (86), in the case of n — s, we find 

3^ 



, T(n + / + 3/2) / 3 
tn = ^(n + Z + - 



When n < s, then 



iL = _r {n+ i + m£(s-P-iy- (n<sh 



si 



p=0 



(n — p)\ 



while in the opposite case, 
T(n + / + 3/2) 



ns 



S\ 



1>: 



•-1 (s-p-i)! — 



(n > s) , 



where 



We also have 



^ = ^(p + 1) - i>(p + s + 1) - i> (p + s + I + - 
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, I> + / + l/2) , r(n + / + 3/2) . T(l + 3/2) 

Vn-i" ( n -i)! ' n! ' 0s s! 

In this way, we can calculate any coefficient 

C nl = lim \ C nl (v) , 

in the sequence (85). In particular, for n = 0, we get 

r(/ + 3/2) ^ (p-1)! 



I + 3/2 ^ p 2 r(p + / + 3/2) 



Appendix C 

The integral (102) can be presented as the limit 

j, = lim JJM - 4,(Q) 

involving the integral 

J» = I™ t^e-t L™(t) L^\t) dt. 
Jo 

Substituting in the latter the Laguerre polynomials, we come to the expressions similar 
to those in Appendix A. In this way, 

, I> + 2/ + 2) " (_i)p r(p + 2/ + 3 + z/) 
J ^ U) = n\ § pi (n- P )!I> + 2/ + 2) " 1 " "> ■ 

where (p n (x) is the same function as in Appendix A. Also we find 

Taking the limit z/ — > 0, we use the asymptotic equalities 

cp n (-p - 1 - v) ~ (-l) p (n - p - 2)! (p + 1)! i/ 
for < p < n — 2, and 

V9„(-p - 1 - !/) ~ (-l) n (p + 1)! {1 + ty>(p + 2) - ^(p - n + 2)] u} 
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for the case of p — n — 1, n. Here the function ip(x) is defined in Appendix B. With these 
properties we obtain the expression 

f nn = ( n + 2l + 1 V- [2n+ i + 2 (n + / + l)^(n + 2/ + 2)] 

defining the coefficients D nt in Eq. (103), 

7-) _ u - J nn 

nl 2(n + l + l)(n + 2l + l)\ " 
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Table captions 



Table I. The accuracy of different approximations for the ground-state energy e(0, 0), 
for various povers v of the potential: S\ and 62 are the percentage errors of the renor- 
malized expressions in Eq. (79); e* 2 is the error of the self-similar approximant (80); and 
£\vkb is the error of the quasiclassical approximation (73). 

Table II. The accuracy of approximations for several energy levels e(0, /) of the log- 
arithmic potential: S\ and 62 are the percentage errors of the renormalized expressions 
(89); el is the error of the self-similar approximant (90); and e WKB is the error of the 
quasiclassical approximation (83). 
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Table I 



V 


e(0,0) 


ei(%) 


e 2 (%) 


e*(%) 


£\vkb(%) 


0.15 


1.2653 


0.26 


0.05 


-0.49 


-0.73 


0.50 


1.5961 


0.44 


0.07 


-0.02 


-1.2 


0.75 


1.7450 


0.39 


0.06 


0.08 


-1.0 


1.5 


2.0121 


0.08 


0.01 


0.04 


-0.26 


2 


2.1213 














3 


2.2765 


0.33 


-0.05 


0.17 


-0.27 


4 


2.3936 


1.3 


-0.43 


0.63 


-1.3 


5 


2.4924 


2.6 


-1.6 


1.2 


-2.6 


6 


2.5797 


4.3 


-4.0 


1.6 


-4.1 


8 


2.7315 


8.5 


-16 


1.0 


-7.1 


10 


2.8616 


13 


-50 


-4.3 


-9.9 



Table II 



I 


e(0,Z) 


£i(%) 


e 2 (%) 


e* 2 (%) 


£WKB\%) 





0.52215 


2.2 


0.48 


0.37 


-6.4 


1 


0.82150 


0.81 


0.10 


-0.36 


-9.4 


2 


1.0075 


0.47 


0.05 


-0.38 


-9.4 


3 


1.1430 


0.31 


0.01 


-0.36 


-9.2 


4 


1.2495 


0.22 


-0.002 


-0.33 


-8.9 
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